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Abstract 

The Greens library is presented which provides a set of C++ procedures for the 
computation of the (radial) Coulomb wave and Green's functions. Both, the nonrelativistic 
as well as relativistic representations of these functions are supported by the library. 
However, while the wave functions are implemented for all, the bound and free-electron 
states, the Green's functions are provided only for bound-state energies (E < 0). Apart 
from the Coulomb functions, moreover, the implementation of several special functions, 
such as the Kummer and Whittaker functions of the first and second kind, as well as a few 
utility procedures may help the user with the set-up and evaluation of matrix elements. 
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PROGRAM SUMMARY 



Title of program: Greens 
Catalogue number: To be assigned. 

Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland. 
Users may obtain the program also by down-loading a tar-file from a home page at the 
University of Kassel 

http : / / www . physik . uni-kassel . de/~kovalp/ software/greens 

Licensing provisions: None. 

Computer for which the program is designed and has been tested: PC Pentium III, PC Athlon 
Installation: University of Kassel (Germany). 

Operating systems: Linux 6.1+, SuSe Linux 7.3, SuSe Linux 8.0, Windows 98. 
Program language used: C++. 

Memory required to execute with typical data: 300 kB. 

No. of bits in a word: All real variables are of type double (i.e. 8 bytes long). 

Distribution format: Compressed tar file. 

CPC Program Library Subprograms required: None. 

Keywords: confluent hypergeometric function, Coulomb-Green's function, hydrogenic wave 
function, Kummer function, nonrelativistic, relativistic, two-photon ionization cross section, 
Whittaker function. 

Nature of the physical problem: 

In order to describe and understand the behaviour of hydrogen-like ions, one often needs the 
Coulomb wave and Green's functions for the evaluation of matrix elements. But although these 
functions have been known analytically for a long time and within different representations 
[1,2], not so many implementations exist and allow for a simple access to these functions. In 
practise, moreover, the application of the Coulomb functions is sometimes hampered due to 
numerical instabilities. 

Method of solution: 

The radial components of the Coulomb wave and Green's functions are implemented in po- 
sition space, following the representation of Swainson and Drake [2]. For the computation 
of these functions, however, use is made of Kummer's functions of the first and second kind 
[3] which were implemented for a wide range of arguments. In addition, in order to support 
the integration over the Coulomb functions, an adaptive Gauss-Legendre quadrature has also 
been implemented within one and two dimensions. 
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Restrictions onto the complexity of the problem: 

As known for the hydrogen atom, the Coulomb wave and Green's functions exhibit a rapid os- 
cillation in their radial structure if either the principal quantum number or the (free-electron) 
energy increase. In the implementation of these wave functions, therefore, the bound-state 
functions have been tested properly only up to the principal quantum number n ~ 20, while 
the free-electron waves were tested for the angular momentum quantum numbers k < 7 and 
for all energies in the range . . . 10|-Ei s |. In the computation of the two-photon ionization 
cross sections ai, moreover, only the long-wavelength approximation (e ,k ' r ~ 1) is considered 
both, within the nonrelativistic and relativistic framework. 

Unusual features of the program: 

Access to the wave and Green's functions is given simply by means of the Greens library 
which provides a set of C++ procedures. Apart from these Coulomb functions, however, 
Greens also supports the computation of several special functions from mathematical physics 
(see section 2.4) as well as of two-photon ionization cross sections in long-wavelength approxi- 
mation, i.e. for a very first application of the atomic Green's functions. Moreover, to facilitate 
the integration over the radial functions, an adaptive Gauss-Legendre quadrature has been 
also incorporated into the Greens library. 

Typical running time: Time requirements critically depends on the quantum numbers and 
energies of the functions as well as on the requested accuracy in the case of a numerical 
integration. One value of the relativistic two-photon ionization cross section takes less or 
about one minute on a Pentium III 550 MHz processor. 
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LONG WRITE-UP 



1 Introduction 

From the early days of quantum mechanics on, the 'hydrogen atom' has served not only as 
a well-known textbook problem but also as one of the fundamental models in the physics of 
atoms, molecules, and nuclei. When combined with the (atomic) shell model, namely, the 
— analytic — solutions for the hydrogen-like ions help understand most atomic processes in 
Nature, at least qualitatively. For this reason also, the 'hydrogen atom' has found its way into 
quite different fields of physics including, for example, astro- and plasma physics, quantum 
optics or even the search for more efficient x-ray lasers schemes. 

Despite of the success of the hydrogen model, however, the Coulomb problem is not always 
that simple to deal with, in particular, if a relativistic treatment is required. Therefore, 
various program tools have been developed over the years to help with either the analytic 
or numerical manipulation of the Coulomb functions and their matrix elements. For the 
nonrelativistic Coulomb problem, for example, the bound-electron states can be obtained 
from the codes of Noble and Thompson [1], who applied a continued fraction representation 
of the Whittaker functions, Bell and Scott [2], or simply by using the Gnu Scientific Library 
[3]. These functions are incorporated also into a recent library by Madsen and coworkers 
[3], which has been designed to support the computation of the multipole matrix elements 
for circular and linear polarized light. - - Less attention, in contrast, has been paid to the 
relativistic wave functions for which a OPC program is provided only by Salvat et al. [5J. This 
program help integrate the radial equation for any spherical-symmetric potential for both, 
the (one-particle) Schrodinger and Dirac equations and also provides separate procedures to 
compute the Coulomb wave functions. 

Apart from the bound and free-electron wave functions, however, the Coulomb Green's func- 
tions play a similar important role, in particular, if the interaction of atoms with external 
fields is to be studied. In second- and higher-order perturbation theory, for instance, these 
functions help to carry out the summation over the complete spectrum in a rather efficient 
way. But although different analytic representations are known for the Green's functions 
[SI El El ID] , until today, there are almost no reliable codes freely available. 

Therefore, to facilitate the further application of the 'hydrogen atom' in different contexts, here 
we present the Greens library which provides a set of C++ procedures for the computation of 
the Coulomb wave and Green's funtions. In Greens, these hydrogenic functions are supported 
both, within a nonrelativistic as well as relativistic framework. Beside of the various routines 
for the computation (of the radial parts) of these functions, however, we also supply the user 
with a Gauss-Legendre quadrature and a set of special functions to simplify the evaluation of 
matrix elements. — But before we shall present details about the organization of the Greens 
library, in the following section, we first compile the basic formulas from the theory of the 
'hydrogen atom' with emphasize especially to those expressions, which have been implemented 
explicitly. In section 3, later, the program structure will be discussed and how the library is 
to be distributed. This section also lists all user-relevant commands, although not much is 
said here about the underlying algorithms. In most cases, we followed the expressions from 
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sections 2 but care has been taken in order to provide a reliable code for a rather wide range 
of parameters which, sometimes, required quite additional effort. In section 4, we explain how 
(easily) the hydrogenic wave and Green's functions can be accessed not only for a particular 
set of arguments but also for the computation of matrix elements. These examples may serve, 
therefore, also as a test bed for the installation of the code. Section 5, finally, gives a brief 
summary and an outlook into our future work. 

2 Theoretical Background 

Since the theory of the 'hydrogen atom' has been presented at quite many places before (see, 
for instance, the texts of Messiah [11] and Drake p2]), we shall restrict ourselves to rather a 
short compilation of formulas, just enough in order to provide the basic notations and those 
expressions which are implemented in the code. In the next two subsections, therefore, we first 
recall the (analytic) form of the Coulomb wave and Green's functions while, in subsection 2.3, 
these functions are applied to calculate the two-photon ionizations cross sections for linear and 
circular polarized light. In all these subsections, the nonrelativistic and relativistic formulas 
are always presented in turn of each other in order to display the similarities but also the 
differences in the numerical treatment of these functions. Subsection 2.4, moreover, provides 
reference to a few special functions from mathematical physics, which frequently occur in the 
computations of the Coulomb wave or Green's functions and, hence, need to be part of the 
Greens library. 

2.1 Coulomb wave functions 
2.1.1 Nonrelativistic wave function 

In a time-independent external field, the motion of a particle is described by the stationary 
Schrodinger equation 



which, obviously, is an eigenvalue equation for the total energy E of the particle. As known 
from the nonrelativistic Schrodinger theory, the Hamiltonian H just includes the kinetic and 
potential energy of the particle and, thus, takes the farnfl 



S(r) = - 




in the case of a (pure) Coulomb field of a nucleus with charge Z. For such a spherical- 
symmetric potential, of course, Eq. ([!]) and the wave functions tp(r) can be separated 



into a radial and an angular part where, in most practical computations, the angular structure 
of the wave functions is often treated by means of the techniques from Racah's algebra |14j . 
In expression ([3|), n and / denote the principal and orbital angular momentum quantum 
numbers, respectively, while m describes the projection of the z— component of the orbital 

1 Here and in the following, we use atomic units (m e = % = e 2 /4ire = 1) if not stated otherwise. 








(1) 




(3) 



r 
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angular momentum onto the quantization axis and is called the magnetic quantum number. 
The radial part of the wave function, P n i (r)/r, is a solution of the radial Schrodinger equation 



1 d ( 2 d \ 1(1 + 1) 2Z 

1 r 2 — ) - V 2 ; + — + 2E 



dr \ dr 



r 



Ml = o (4) 



which has (normalizable) physical solutions for a discrete set of negative energies 

Z 2 

E n = < 0, n = 1, 2, ... , (5) 

the so-called bound states, as well as for all positive energies E > 0, i.e. the continuum or 
free-electron states. Both, the bound and continuum solutions of (jH) can be represented in 
terms of a single Whittaker function of the first kind M a ^ {z) 

P nl (r) = C(n,l,Z)M n , + 1 (2Zr/n) (6) 



P El (r) = C(E,l,Z)M. ,(-2ix/2#Ir) (7) 



iJ^- 1+ 



with real or complex arguments, and where C(n,l,Z) and C(E,l, Z), respectively, denote 
the corresponding normalization factors. The Whittaker functions are closely related to the 
Kummer functions of first and second kind as we will discuss in subsection 2.2.1. In the 
standard theory, moreover, the radial wave functions © and ([7]) are often normalized due to 



oo 

2 



Pnl(r)dr = 1, (8) 
P* El {r)P m {r)dr = 5(E — E') , (9) 



o 



in order to represent a single particle per bound state or per energy unit, respectively, if 
particles in the continuum are concerned. 

2.1.2 Relativist ic wave functions 

An eigenvalue equation analogue to dU also applies, if the motion of the particle is described 
within the relativistic theory. For an electron with spin s = 1/2, however, then the Hamilto- 
nian H needs to be replaced by the Dirac-Hamiltonian [11] 

H D ( r ) = -ica • V + 13c 2 - — (10) 

r 

which, apart from the kinetic and potential energy of the electron in the field of the nucleus, 
now also incorporates the rest energy of the electron as well as energy contributions owing to 
its spin. As in the nonrelativistic case, a separation of the wave function 

, , , I ( Pn^{r) tt Km (0,<p) \ ni , 
^/W(r) = - (11) 

into a radial and angular part is possible for any spherical-symmetric potential, where the two 
radial functions P nK (r) and Q nK (r) are often called the large and small components. These 
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two functions also form a radial spinor 
first-order, coupled equations [13] 



Pnn (r) 

QnK iy , 



and have to be obtained as solutions of the 



-— - E 

r 



Id k 

— H 

ar or ar 



Pnn jr) 

r 

Pnn (r) 

r 



+ 



K 

ar 



l_d_ 

ar dr 



2 z 
-77 + - + E 



a*- 



QnK (jj 



QnK (r) 







(12) 



(13) 



where, however, the (total) energy E is taken here to represent the energy of the electron 
without its rest energ>@ c 2 , similar to Eq. ^ in the Schrodinger theory. In Eqs. (|llH13p . 
moreover, k = ±(j + 1/2) for I = j ± 1/2 is called the relativistic angular momentum 
quantum number and carries information about both, the total angular momentum j as well 
as the parity (— 1)' of the wave function. Again, (normalizable) physical solutions to the Dirac 
operator (jlOp can be found for a discrete set of negative energies 



E n 



a 



1 + 



aZ 



n 



a 



2 Z 2 



n 



1, 2, 



-n, 



a -2 < 



n 



1, K ^ 



(14) 



and for all positive energies E > as well as for the (negative) energies E < —2c 2 . The 
two latter — contineous — parts of the spectrum are also called the positive and negative 
continuum whereby the negative branch, in particular, requires some re-intepretation of the 
theory (in terms of positron states, for example) and often introduces additional complications 
in the treatment of many-electron systems. When compared with the nonrelativistic energies 
([5]), however, the degeneracy of the (relativistic) energies (|14p is partially resolved and now 
depends on both, the principal quantum number n and the relativistic quantum number k. 

Explicit representation of the bound and free-electron solutions of Eqs. (|12H13|) are known 
from the literature (cf. [13, [15]) but typically result in rather lengthy expressions. For the 
bound states, for example, the two radial components are given by 



Pn 



QnK (jj 



C P (n,K,Z) r s e~ qr [(-n + \k\) M(-n + |«| + 1, 2s + 1; 2qr) 

- (k - Zq- 1 ) M(-n + \k\, 2s + 1; 2qr)] 

C Q (n,K,Z) r s e~ qr [-(-n + \k\) M(-n + \k\ + 1, 2s + 1; 2qr) 

- (k -Zq' 1 ) M(-n + 2s + 1; 2qr)) 



(15) 



(16) 



where M(a, b; z) is the Kummer function of the first kind, s = \/ n 2 — (aZ) 2 , and q = 
Z [(aZ) 2 + (n — | At | +s 2 )] - 2 ; while even more elaborate expressions arise for the free-electron 
states [15]. Similiar to ([8]) and ([9|),the bound and free-electron radial wave functions can be 
normalized also due to 

f oo 

{Pl(r) + Ql(r)) dr = 1, (17) 



POO 

/ (Pek {r)P E > K (r) + Qek {t)Qe'k (r)) dr = 5{E - E') 
Jo 



(18) 



In atomic units, the speed of light c = 1/a is the inverse of the fine-structure constant. 
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to represent one electron per bound state or per energy unit, respectively. 

In the Greens library, the radial functions of the bound and free-electron states can be ac- 
cessed by means of the two library procedures greens_radial_orbital() and greens_radial_spinor() 
in the nonrelativistic and relativistic case, respectively; for further details, see section 3. 

2.2 Coulomb Green's functions 

Apart from the wave functions, which describe the electron in particular quantum states, 
one often needs a summation over all (unoccupied) states, especially, if parts of the atomic 
interaction are treated as a perturbation. A full summation is required in second- and higher- 
order perturbation theory, for instance, if the behaviour of the atom is studied in a — not too 
weak — radiation field or in the presence of external electric or magnetic fields. Although, 
in principle, it appears straightforward to carry out such a summation explicitly, the large 
number of terms and the need of free-free matrix elements may hamper such an approach. 
Instead, the use of Green's functions [16] 

G£(ry) = % i*w>y>i m 

often provides a much simpler access to the spectrum of the atom and, hence, to a perturbative 
treatment of atomic processes. In the following, therefore, we first recall a representation 
of the radial Coulomb Green's functions as appropriate for numerical computations. The 
application of these functions in the computation of two-photon ionization cross sections 02 
for hydrogen-like is discussed later in subsection 12.31 



2.2.1 Nonrelativistic Green's function 

Analogue to the wave functions ([3]), the Coulomb Green's functions G^(r, r') are obtained as 
solutions of a linear equation 

(H(r)-E)G E (r,r') = 6(r - r') (20) 

with the same Schrodinger operator as in (pQ) but for an additional 5 — like inhomogenity on 
the right-hand side, which allows for solutions for any arbitrary E. For a spherical-symmetric 
potential, again, this equation can be separated into a radial and angular part by using the 
ansatz 

G E (r,r') = Y, gEl YU^)YL(AV) (21) 

lm 

for the Green's function in spherical coordinates. By substituting ansatz (I2ip into Eq. (|20l) . 
one easily shows that the radial Green's function g E i (r, r'), which just depends on the energy 
E and the orbital angular momentum /, must satisfy the equation 



1 d ( 2 d \ 1(1 + 1) 2Z 

1 r 2 — - K „ ; + — + 2E 



Iei (r, r') _ _ 2 5{r - r') 



Solutions to this single equation can be determined by taking a proper superposition of the 
regular and irregular solutions (near the origin) of Schrodinger 's equation @. An explicit 
representation for the radial Green's function reads as [9] 

gra(r ' r,) = T xr^a + 2) M ^ (2xr<) w ^ (2xr>) ' (22) 
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maxfr, r') and r . 



minfr, r') refer to the 



where x = (— 2E) I , r = — , and where r> - iuhaw./ / aim 

x 

larger and smaller value of the two radial coordinates, respectively. In this representation, 
moreover, M a ^ (z) and W aj 6 (z) denote the two Whittaker functions of the first and second 
kind which can be expressed also in terms of the Kummer functions M(a, b; z) and U(a, b; z) 
of the corresponding kinds [T7] 

M a , 6 (z) = z b+ ^e- z/2 M(b-a + 1/2, 2b + l;z), (23) 

W a , 6 (z) = z b+ ^e- z/2 V(b-a + 1/2, 2b + 1; z) . (24) 

In practise, the two Kummer functions are used more frequently (than the Whittaker func- 
tions) in the mathematical literature and in various program libraries since M(a, b; z) is closely 
related to the hypergeometric series and since the Kummer function U{a,b;z) of the second 
kind can be expressed in terms of M{a,b;z). In addition, several improved algorithms have 
been worked out recently in order to calculate the regular Kummer function M{a, b; z) more 
efficiently, see section [2~4l for further details. 

2.2.2 Relativistic Green's function 

Of course, the relativistic Coulomb Green's function must refer to the Dirac Hamiltonian (jlOp 
and, hence, is given by a 4 x 4-matrix which satisfies the equation 



H D (r)-E-c 2 ) G E (r,r' 



= <5(r-r')l4, 

where 1 4 denotes the 4x4 unit-matrix and where, as for the wave functions from Eqs. (|12| - 
[T3j) . the rest energy c 2 has not been incorporated into the (total) energy E. Solutions to 
this equation are known again from the literature for a radial- angular representation of the 
Coulomb Green's function [9] 



_ 1 / ^(r,r')l] Km (r)!]t m (r') 
G E r,r') = E — 

rr 1 ^gI K L (r,r')^ Km (r)^^t m (r') 



-tg|f (r, r')Q Km (v)Ql Km (v') 
(r, r') n. Km (r) nl Km (r<) 



r ss 



(25) 



, , /g|*(ry) g|f(r,r') 
where the radial part 

which must satisfy the matrix equation 



/rr' of this function is now a 2 x 2-matrix 



/ 



[-— - E 




r 





1 d k 

— H 

ar or ar 



K 

ar 



OL 



J_d_ 

ar dr 

Z 

r 



E 



\ 



'g|«(r,r') g^(r,r'Y 

k &Ek ( r > r ') g|f ( r > r ')j 



5{r 



In this representation of the Green's function, we make use of the two superscripts T and T' 
to denote the individual components in the 2x2 radial Green's matrix. They may take both 
the values T = [L, S] to refer to either the large or small components, when multiplied with 
a wave function spinor (jlip . An explicit representation of the (four) components g^J (r, r') 
of the radial Green's function is found by Swainson and Drake [9] 



' h 11 - X(h 12 + h 21 ) + X 2 h 22 
-X(h n + h 22 ) + X 2 h 12 + h 21 



-X(h n + h 22 ) + h 12 + X 2 h 2ls 
X 2 h n - X(h 12 + h 21 ) + h 22 , 



(26) 



9 



with 



«"M - (l^!ffi^±W - ') r(7->) M„_,(^)W„ T _,(2^) ,28) 



fc2i(r, r') = /i 12 (r', r) 

(1 - X 2 )T{ 1 + 1 - i/) cry" 1 
2T(2 7 + 2) 



2 7 (2 7 + 1) 0(r' - r) M !/j 7 _i (2wr) W w , 7+ i (2a;/) 
- (z/ + 7 )e(r-/)W V!7 _i(2 W r)M !/)7+ i(2a;r / )l (29) 



and 



X = (-K + 7 )(aZ)- 1 , 7 = (k 2 - a 2 ^ 2 ) 1 / 2 , 

to = cT l (l- {Ea 2 + l) 2 ) l/2 , u = Z(Ea 2 + 1)uj- 1 , 
and where 9{x) denotes the Heaviside function. 

In the Greens library, we provide the two procedures greens_radial_function() and 
greens_radial_matrix() which support the computation of the radial functions (|22p and ([26 
for any proper set of parameters. 

2.3 Two-photon transition amplitudes and ionization cross sections 



The Green's function (|2ip and (J25J) can be utilized directly to evaluate, for instance, the two- 
photon cross sections o<i for a non-resonant excitation, ionization, or decay process. They 
also occur rather naturally in the theory of the photon scattering on hydrogen-like ions. In 
the following, we briefly outline the perturbative calculation of the two-photon ionization 
cross section for hydrogen-like ions which, for an unpolarized target and in atomic unitjf], is 
given by 



8vr 3 a 2 

7 Ktm 



E^ViEi^ (30) 



£ 2 ^ 2j i + 



where E~ is the photon energy and Mfi the two-photon transition amplitude 

■ ik 2 r . r> I ih..\ Ilk. I 1^ P ik l r 



M 



% 4rf E v — E^ — Ei 

In this amplitude, moreover, (tpi, Ei), (tp u , E u ), and (tpf, Ef) denote the wave functions and 
energies of the initial, intermediate and final atomic states, respectively. Here, the energy of 
the final state, Ef, does not appear explicitly in (f3Tj) but follows from 

Ej = E{ + 2E~/ 

due to the conservation of energy. Furthermore, the two vector quantities u\ and p in the 
transition amplitude (f!?Tj) refer to the polarization of the two photons as well as to the electron 
momentum operator. 



3 The cross section <J2 has the dimension length 4 x time and, thus, can be converted into cgs-units cm 4 -s by 
using the multiplication factor 1.896792 ■ 10~ 50 . 
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As mentioned before, the summation over v in (|3ip runs over the complete spectrum of the 
atom including the continuum. This summation can be replaced, therefore, by a single Green's 
function (fT9|) . so that the transition amplitude ([3T|) finally takes the form 

M fi = y"^(r)u A2 e ik2r -pG^ +i?7 (r,r')u Al e ikir '-p'^(r , )rfr(ir'. (32) 

It is this form of the transition amplitude which has often been used in the literature to study 
non-resonant, two-photon processes [TU1 [T8]. 

2.3.1 Nonrelativistic ionization cross sections 

For the sake of brevity, let us restrict ourselves to the two-photon ionization cross sections 
within the long-wavelength approximation, i.e. we assume e lkr = 1 for the coupling of the 
radiation field in (|32|) . Apart from the electric-dipole field, of course, this approximation 
neglects the contribution from all higher multipoles, but is known to describe well the ioniza- 
tion of light atoms with a nuclear charge of, say, Z < 30 and for photon energies below the 
ionization threshold E~ < Et- By substituting p —* r and E^ —* into Eqs. (|30p and 

(|31|) . moreover, we may obtain the ionization cross section in length gauge 



8 - 3 «X E ^£i< ngth) i 2 > ( 33 ) 



(length) 



l f m f 



with 



J V} (r) u A2 • r G El+El (r, r') u Al • r' ^ (r') dr dr' . (34) 



M (length) 

Using the radial- angular representations ([3]) and (|2ip of the wave and Green's functions, re- 
spectively, and by making use of some angular momentum algebra, the 6-dimensional integral 
in the transition amplitude (|34p can be reduced further to just a two-dimensional integration 
over the radial coordinates r and r'. In addition, if we assume the ion initially in its Is ground- 
state and circular polarized light, i.e. two photons with the same helicity Ai = A2 = ± 1, the 
two-photon ionization cross section (in length gauge) simply takes the form 



^length, circular) = ^2 ^2 



J Pe;2 (r) r g£ l+ £ 7 ,i (r, r') r' P 10 (/) dr dr' 



(35) 



2.3.2 Relativistic two photon ionization cross sections 



The long-wavelength approximation for the coupling of the radiation field can be considered 
also within the framework of the relativistic theory. In this framework, however, an useful 
estimate of the total cross section 02 are obtained only if the photon energy is well below 
the threshold energy E~ < Et of the two-photon ionization. In the relativistic theory, the 
(long-wavelength) transition amplitude (|32|) takes the form 

M fi = c 2 J ^/(r)u A2 -cxG^+e, (r,r')u Al ■ a! ^(v 1 ) dv dr' , (36) 

where a denotes Dirac's velocity operator. Using the radial-angular representation (|25|) of the 
Green's functions, then the total two-photon ionization cross section o~ 2 for circular-polarized 
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light can be written as 



(velocity, circular) 87T ( ->2 



3C/ 5L (d3, P 3, S i) - 5[/ Li (d3,pi, S i) - 15[/ i5 (d3,p., s .)] I, (37) 



where we introduced the radial integral 

U TT '( Kf , k u , m) = J gl fKf (r) g^ +EytKv (r, r') (r') dr dr' . (38) 

In this integral, a superscript T refers to the conjugate of T, i.e. T = S for T = L and 
wee versa, and gn K ( r ) an d gn K ( r ) are used to denote the large and small components of 
the radial spinor pip. This notation allows for a very compact representation of the multi- 
photon transition amplitudes which can be applied also well beyond the long-wavelength 
approximation . 

In the Greens library, the procedure greens_two_photon_cs() is presented to compute two- 
photon ionization cross sections in various approximations. 



2.4 Special functions 

Of course, the main emphasize in developing the Greens library has been paid to the compu- 
tation of the Coulomb wave and Green's functions as appropriate for a theoretical description 
of hydrogen-like ions. As seen from sections 2.1 and 2.2, however, for an explicit represen- 
tation of these functions we usually need to refer to a few special functions such as the T(z) 
and ^f(z) functions, or the Kummer and Whittaker functions of the first and second kind 
which are known from the mathematical literature [17] ■ Therefore, in order to facilitate the 
implementation of the Coulomb functions, we have to provide also a simple interface to these 
special functions; in the following, we briefly summarize the definition of these functions and 
for which type of arguments they are needed for the Greens library. 

Euler's Gamma function T(z) and the Psi-function ^f(z) occur very frequently and in quite 
different fields of physics. While the T— function is defined by the integral 

/>oo 

F(z) = / f^e-'alt (39) 
J o 

the ^—function refers to the derivative 

V(z) = 1 ^ n . (40) 

These functions are defined for all complex arguments z except of the real negative integers z ^ 
— 1, —2, ... where they have their poles. In Greens, the T(z) function with real arguments 
z is needed for the computation of the bound-state wave and Green's functions, respectively, 
while complex arguments arise in the representation of the free-electron waves (J7J)- The 
^—function, in addition, arises in the calculation of the Kummer function U(a, b; z) of the 
second type if the argument b refers to an integer in the computation of nonrelativistic Green's 
functions. 
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Although the Coulomb wave and Green's functions are often expressed in terms of the Whit- 
taker functions M aj 6 (z) and W a j (z) of the first and second kind, in practical computations 
one makes better use of the Kummer functions of the corresponding kind, as discussed in sub- 
section 2.2.1 above. The Kummer functions M(a,b;z) and U(a,b;z) of the first and second 
kind refer to the regular and irregular solutions of Kummer's equation 



d 2 M „ . dM 

z ^ + (i, - 2) d7- aM = 0; (41) 



in the literature, however, also several other notations are used for these functions such as 
Mfa, b; z) = \Fi(a; b; z) or Ufa, b, z) = ^(a, b, z), respectively. Usually, the function M(a, b; z) 
of the first kind is solved for the initial value M(o, b; 0) = 1 and, hence, is given by the confluent 
hypergeometric series 

, w , % o 1 of a + 1) 9 . , . 

M(„,6;,) = l + ^ + -^L TI i^ + .... (42) 

The Kummer function of the first kind Mfa, b; z) is needed for both, real a, b, z and complex 
arguments a, z to represent the radial wave and Green's function components. In contrast, 
the Kummer function Ufa, b; z) of the second kind is required only for real argument b, for 
which it can be expressed as a linear combination 



7T 

Ufa, b; z) 



sin irb 



M(a,b;z) x _ h M(a + 1 - b, 2 - b; z) 



r(i + o-6)r(6) " r(o)r(2-6) 



(43) 



of two Kummer functions of the first kind; the function Ufa, b; z) arises in the computation 
of the radial Green's function. 

The following section explains how these special functions from the Greens library can be 
used also in applications other than the computation of Coulomb wave and Green's functions. 



3 Program organisation 

3.1 Overview about the Greens library 

The Greens library has been designed mainly in order to facilitate numerical applications 
of the Coulomb wave and Green's functions from section 2. It provides the user with a 
set of C+- 1- procedures to compute the radial components of these functions within both, a 
nonrelativistic as well as relativistic framework. Apart from the radial components, however, 
we also support the numerical integration of the Coulomb functions as well as the computation 
of a few selected matrix elements which, below, will help us demonstrate the application of 
the Greens library. To provide the user with a simple access to the various functions, the 
concepts of object-oriented programming such as structures, classes and members as well as 
the overloading of procedures and operators have been utilized carefully. 

Table 1 lists the main procedures of the Greens library for calculating the energies and 
radial components of the Coulomb functions. To simplify the use of the library, the classes 
spinor2_col, spinor2_raw, and matrix_2x2 have been implemented to describe the radial 
spinor (jlip . its adjunct raw spinor, and the radial Green's matrix (|26p . respectively. The 
classes spinor2_col and spinor2_raw, for instance, contain each the two members .L and 
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.S to represent the large and small components of a relativistic wave function, while the 
class matrix_2x2 has the four members . LL, .LS, .SL, and .SS with an obvious meaning. 
The class matrix_2x2, moreover, also contains the member . e which just returns all the four 
matrix elements together within a 2 x 2 array. 

In order to treat the bound- and free-electron states in a similar way, the two wave func- 
tion procedures green_radial_orbital() and green_radial_spinor() have been overloaded. For 
these two procedures, a first integer argument n > 1 is used to represent the principal 
quantum number and to return the corresponding bound-state solution, while a (first) argu- 
ment E > of type double refers to the kinetic energy of a free-electron state (in Hartree 
units). As mentioned above, however, this energy E does not include the electron rest en- 
ergy, neither in the nonrelativistic nor relativistic framework. The two additional procedures 
greens_set_nuclear_charge() and greens_get_nuclear_charge() from table 1 can be called to re- 
define or to return the current value of the nuclear charge which is utilized for the computation 
of all radial functions. The default value of the nuclear charge is Z = 1. 

In most applications, the (radial) Coulomb wave and Green's function components usually 
occur as part of some matrix element and, hence, first require an additional integration (over 
r and/or r') before any observable quantity is obtained. Therefore, to facilitate such appli- 
cations, we also provide the utility procedure greens_integral_GL() which evaluates a 1- or 
2-dimensional integral over a finite or infinite area with a user-defined accuracy, see table 2. 
In this procedure, a Gauss-Legendre quadrature [T7] of appropriate order is applied, indepen- 
dently for each dimension of the integrand. Moreover, to ensure a result which is accurate 
up to a given number of d valid digits, the domain of integration is divided by steps into 
subdomains until the required accuracy is obtained. A Warning arises during the execution, 
if the requested precision cannot be guarranteed by the procedure. As seen from table 2, the 
procedure name greens_integral_GL() is overloaded and, thus, can be invoked with rather 
different lists of parameters, from which the dimension of the integral, the integration domain 
as well as the type of the function is deduced. Apart from a real-valued integrand fix) or 
f(x,y), respectively, greens_integral_GL() also supports a matrix_2x2-valued integrand as 
appropriate for the computation of matrix elements such as (|38p from the relativistic theory. 
In the latter case, for instance, all the four integrals U LL , U LS , U SL and II s s in (|37p could 
be treated simultaneously. 

A second utility procedure greens_two_photon_cs() from table 2 enables the user to calculate 
two-photon ionization cross sections in various approximations. Obviously, this procedure 
makes use of greens_integral_GL() and is mainly provided for test purposes below. It helps 
compute the total two-photon ionization cross sections 02 for linear or circular polarized light 
and within either the nonrelativistic or relativistic framework, respectively. In all of these 
cases, however, the computation of the cross sections is restricted to the long-wavelength 
approximation e jkr = 1 for the coupling of the radiation field and to the ionization of an 
electron from the unpolarized Is ground state. In addition, the photon energy E~, i.e. the 
third argument of the procedure greens_two_photon_cs() must be in the range —E\ B /2 < 
E-y < —Eis where Ei s is the (negative) Is— binding energy from Eqs. © or ()14|) . Again, the 
last argument d refers to the requested accuracy of the cross section of (at least) d valid digits 
and is transfered directly to the underlying integration procedure greens_integral_GL() . 

Of course, the wave and Green's functions from section 2 can hardly be implemented without 
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Table 1 Main procedures of the Greens library to calculate the energies and radial wave and Green's functions for hydrogen-like ions. The (expected) 
type of parameters is shown by using the syntax of CH — h; all quantities below must be given in atomic units. 



Procedure 



Arguments 



Description and comments 



double greens_energy (int n) 

(int n, int kappa) 

double (int n, int 1, double r) 

greens jradial_orbital 

(double E, int 1, double r) 

spinor2_col (int n, int kappa, double r) 

greens jradial_spinor 

(double E, int kappa, double r) 



double (double E, int 1, double r, 

greens jradial Junction double r') 

matrix_2x2 (double E, int kappa, double r, 

greens jradialjnatrix double r') 



Returns the nonrelativistic energy E n (in a. u.) of a bound-state solution with 
principal quantum number n; see Eq. ([5]). 

Returns the relativistic energy E nK (in a. u.) of a bound-state solution quantum 
numbers n and k; see Eq. (fH)) . 

Computes the value of the radial function P n i (r) at r of a bound state ^ with 
principal quantum number n and orbital angular momentum I. 
Computes the value of the radial function Pei (r) at r of a free-electron state ([7]) 
with energy E > and orbital angular momentum I. 

fPn K (r)\ 

Computes the value of the radial spinor function at r of a bound state 

\Q nK (r)J 

(|15H16|) with principal quantum number n and relativistic angular momentum 
quantum number k. 

Computes the value of the radial spinor function at r of a free- 

\Q EK .(r)J 

electron state with energy E and relativistic angular momentum quantum number 

K. 

Computes the radial Coulomb Green's function gEi (r, r') at r and r' from ([22]) 
for the energy E < and orbital angular momentum I. 

(g g^"^ \ 
co from (|26p at r and 
g|« SekJ 

r' for the energy E < and the relativistic angular momentum quantum number 
K. 



Table 2 Utility procedures of the Greens library for the numerical integration of 1- and 2-dimensional functions and the computation of two-photon 
ionization cross sections a 2 hi various approximations. The same notation as in table 1 is used. 



Procedure 



Arguments 



Description and comments 



double 

greens_integral_GL 



double 

greens_two_photon_cs 



(double (*funct) (double x) , double a, 
double b , int d) 



(double (*funct) (double x) , int d) 



(double (*funct) (double x, double y) , 
double ax, double bx, double ay, 
double by, int d) 

(double (*funct) (double x, double y) , 
int d) 



( "nonrelativistic" , "circular" , 
double E_ph, int d) 



( "nonrelativstic" , "linear", 
double E_ph, int d) 
("relativistic" , "circular", 
double E_ph, int d) 

("relativistic", "linear", 
double E_ph, int d) 



Calculates the definite (1-dimensional) integral J f(x) dx with an accuracy of 

a 

(at least) d valid digits. This procedure applies an adaptive Gauss-Legendre 

integration formula, independently in each dimension. 

00 

Calculates the definite (1-dimensional) integral J f(x)dx with an accuracy of 



(at least) d valid digits if f(x) does not oscillate rapidly and vanishes sufficiently 
fast for large values of x. 

Calculates the definite (2-dimensional) integral J J f{x, y) dxdy with an accu- 
racy of (at least) d valid digits. 

00 00 

Calculates the definite (2-dimensional) integral J J f(x,y) dxdy with an accu- 



racy of (at least) d valid digits if f(x,y) does not oscillate rapidly and vanishes 
sufficiently fast for large values of x and y. 



Computes the nonrelativistic two-photon ionization cross section Q35[) for circu- 
lar polarized light, in long- wavelength approximation, and for a photon energy 
-Eph > Eis/ '2. A cross section value in atomic units and with an accuracy of (at 
least) d valid digits is returned. 

Computes the nonrelativistic two-photon ionization cross section for linear po- 
larized light and in long-wavelength approximation. 

Computes the relativistic two-photon ionization cross section (|37[) for circu- 
lar polarized light, in long- wavelength approximation, and for a photon energy 
E p h > Eis/ '2. 

Computes the relativistic two-photon ionization cross section for linear polarized 
light and in long-wavelength approximation. 



Table 3 Special function procedures of the Greens library. The same notation as in table 1 is used. 
The type of all procedures is double if all arguments are double, and is of type complex otherwise. 

Procedure Arguments Description and comments 

GAMMA (double z) or (complex z) Returns the T(z) function (f39|) for either a 

real or complex argument z. 

Psi (double z) or (complex z) Returns the ^f(z) function (jlQ)) for either a 

real or complex argument z. 

KummerM (double a, double b, double z) or Calculates the Kummer function M(a, b;z) 

(complex a, double b, complex z) of the first kind (|42|) for real and/or complex 

arguments a, b, and z. 

KummerU (double a, double b, double z) Calculates the Kummer function U(a, b; z) of 

the second kind (pt3"|) for real arguments a, b, 
and z. 

WhittakerM (double a, double b, double z) or Calculates the Whittaker function M aj b(z) of 
(complex a, double b, complex z) the first kind (|23p for either real or complex 

arguments a, 6, and z\ b must be real. 
WhittakerW (double a, double b, double z) Calculates the Whittaker function W a ,b(z) of 

the second kind ([24]) for real arguments a, b, 

and z. 



a proper set of special function procedures. Therefore, table 3 displays those procedures which 
are provided by the Greens library and which we briefly discussed in section 2.4. The allowed 
types of the parameters are also displayed in this table. 

3.2 Distribution and compilation of the Greens library 

The Greens library will be distributed as the gzipped tar-file greens. tar. gz from which 
the greens root directory is obtained by gunzip greens. tar. gz and tar -xvf greens. tar. 
This root contains a Read. me file, the src subdirectory for the source code as well as six 
subdirectories for various examples. In src, we provide the header file greens. h and a 
makefile to facilitate the compilation of the (static) library libgreens.a in the greens root 
directory. It also incorporates about 50 source files for all of the individual procedures. 

In the following section, two examples from the subdirectories example-coulomb-f unct and 
example-twophoton-cs are discussed in more detail and are taken as the test for the instal- 
lation of the library. Each of these example subdirectories, again, contain a makefile from 
which an executable (a. out) is generated simply by typing make within the corresponding 
subdirectory. Since these makefiles also compile and link the library libgreens.a, the user 
may start directly from a copy of one of these subdirectories for his own application of the 
Greens library. 

4 Examples 

To illustrate the use of the Greens library, we first show how the (radial) Coulomb wave 
and Green's functions can be calculated for any point r or (r, r'), respectively. Hereby, 
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#include "greens. h" 



int main (void) { 

int n, 1, kappa; // quantum numbers 

double r, rp, E, wf_nr, gf_nr; // coordinates, energies, etc. 

spinor2_col wf_r; // relativistic spinor 

matrix_2x2 gf_r; // relativistic Green's matrix 

print("#Test of the Coulomb radial functions"); 

for(double Z=1.0; ZO3.0; Z=Z+91.0) { 
print () ; 

greens_set_nuclear_charge (Z) ; // set nuclear charge 

E = -greens_energy (1) * 0.8; 

rp = 2.5/Z; n = 4; 1=2; kappa = -3; 
write("# coord wf_nr wf_r.L "); 

print (" gf_nr gf_r.e[0][0] gf_r.LL"); 

for (r=0.0; r<25.0/Z; r=r+0.1/Z){ 

wf_nr = greens_radial_orbital(n, 1, r) ; 

wf_r = greens_radial_spinor (n, 1, r) ; 

gf_nr = greens_radial_f unction(-E, 1, r, rp) ; 

gf_r = greens_radial_matrix (-E, 1, r, rp) ; 

printf("7„E 7.E 7„E 7.E 7.E 7.E\n" , r, wf_nr, wf_r.L, 

gf_nr, gf_r.e[0] [0] , gf_r.LL); » 

return ; } 

Figure 1: Calculation of the Coulomb wave and Green's functions for nuclear charge Z = 1 
and Z = 92. The printout of this procedure is shown in the Test Run Output and in the 
file printout.txt in the subdirectory example-coulomb-funct. 

a simple comparison between the nonrelativistic and relativistic theory — in the limits of 
a low and high nuclear charge Z — is achieved by setting Z = 1 (hydrogen) and Z = 
92 (hydrogen-like uranium), respectively. Figure 1 displays the source code which evalu- 
ates the two radial functions P/±d{r) and i-4d 5/2 (r), respectively, for r— values in the range 
r = 0., . . . , 25. jZ with a stepsize of Ar = 0.1/Z. Beside of these wave function com- 
ponents, this code also calculates the Coulomb Green's functions at the same values of r 
and for a fixed r' = 2.5/Z. For a call of this procedure, the printout is (partially) shown 
in the Test Run Output below. The source of this example and the complete print- 
out can be found in the subdirectory example-coulomb-funct. In order to obtain the 
— full — radial part of the Coulomb wave and Green's functions, of course, the results of 
greens_radial_orbital() and greens_radial_spinor() must be multiplied with 1/r, while the val- 
ues from greens_radial_function() and greens_radial_matrix() have to be multiplied with 1/r r', 
respectively. 

A second example concerns the computation of the two-photon ionization cross sections for 
the two ions from above. For these ions, the Is binding energies are —1/2 and —4232 Hartrees 
within the nonrelativistic theory. In the Test Run Output below, the two-photon ionization 
cross sections for circular and linear polarized light and within both, the nonrelativistic and 
relativistic approximation. For each of these ions, the cross sections are calculated with an 
accuracy of about six digits for the ten energies E Q , E Q + 0.01 * Z 2 , . . . , E a + 0.09 * Z 2 
where E Q = 0.3 * Z 2 corresponds to 60 % of the nonrelativistic Is binding energy. Again, 
the full source of this example is provided with the Greens library in the subdirectory 
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example-twophoton-cs and, thus, can easily be modified for any other photon energy. 



5 Summary and outlook 

To facilitate applications of the 'hydrogen ion model' in quite different fields of physics, the 
Greens library is presented and provides a set of C++ procedures for the computation of 
the Coulomb wave and Green's functions within both, a nonrelativistic as well as relativistic 
framework. Since C++ is today freely available for most architectures, an object-oriented 
approach to the Coulomb problem could be realized without the need for special compilers 
or other mathematical libraries. Apart from the radial Coulomb functions, however, Greens 
also provides a set of special functions as well as a few utility procedures to evaluate, for 
instance, the two-photon ionization cross sections in long-wavelength approximation. 

In the future, various extensions of the Greens library might be of great interest for the 
physics community. Owing to the current design of several free-electron laser (FEL) facilities 
worldwide, for example, systematic investigations on multiphoton processes become more and 
more likely also in the EUV and x-ray region, where the inner-shell electron get involved. For 
such investigations, which will consider also many-electron atoms and ions, the generation of 
effective one-particle Green's functions are certainly desirable. First steps into this direction, 
including the combination with the well-known Ratip package [19], are currently under work 
in our group. 
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Test Run output 



A. Computation of the radial Coulomb wave and Green's functions 

#Test of the Coulomb radial functions 



#Nuclear charge is changed to 1.000000 
# coord wf_nr wf_r.L 

0.000000E+00 0.000000E+00 0.000000E+00 
1.000000E-01 6.758392E-06 6.761076E-06 9 
2.000000E-01 5.228909E-05 5.230066E-05 7 

2.480000E+01 -2.307241E-01 -2 . 307229E-01 
2.490000E+01 -2.295400E-01 -2 . 295388E-01 



gf_nr gf_r.e[0][0] gf_r.LL 

.000000E+00 0.000000E+00 0.000000E+00 

.861696E-05 9.865164E-05 9.865164E-05 

.642130E-04 7.643481E-04 7.643481E-04 

8.979656E-09 8.980284E-09 8.980284E-09 
8.243999E-09 8.244579E-09 8.244579E-09 



#Nuclear charge is changed to 92.000000 
# coord wf_nr wf_r.L 

0.000000E+00 0.000000E+00 0.000000E+00 
1.086957E-03 6.482422E-05 4.344695E-04 1 
2.173913E-03 5.015393E-04 1.952599E-03 8 

2.706522E-01 -2 . 201670E+00 -2 . 088423E+00 
2.717391E-01 -2.190133E+00 -2 . 074839E+00 



gf_nr gf_r.e[0][0] gf.r.LL 

.000000E+00 0.000000E+00 0.000000E+00 

.071923E-06 5.170738E-06 5.170738E-06 

.306663E-06 2.311893E-05 2.311893E-05 

8.960868E-11 1.641947E-10 1.641947E-10 
8.226637E-11 1.512671E-10 1.512671E-10 



B. Computation of two-photon ionization cross sections 

#Test of the two-photon ionisation cross sections 
#Digits is changed to 6 



#Nuclear charge is changed to 1.000000 

# E cs_nr_c cs_r_c cs_nr_l cs_r_l 

3.000000E-01 8.728681E-01 8.727765E-01 5.849625E-01 5.849002E-01 

3.100000E-01 8.819793E-01 8.818399E-01 5.889291E-01 5.888364E-01 

3.200000E-01 9.143732E-01 9.143980E-01 6.095973E-01 6.096138E-01 

3.800000E-01 5.778480E+00 5.794808E+00 4.829547E+00 4.843563E+00 

3.900000E-01 1.792990E-01 1.796899E-01 2.330477E-01 2.333333E-01 



#Nuclear charge is changed to 92.000000 

# E cs_nr_c cs_r_c cs_nr_l cs_r_l 

2.539200E+03 1.439533E-12 6.927729E-13 9.647196E-13 4.629667E-13 

2.623840E+03 1.454559E-12 6.763950E-13 9.712612E-13 4.510878E-13 

2.708480E+03 1.507983E-12 6.588038E-13 1.005347E-12 4.392479E-13 

3.216320E+03 9.529863E-12 6.744809E-13 7.964884E-12 4.646710E-13 

3.300960E+03 2.956996E-13 7.498654E-13 3.843421E-13 5.232579E-13 
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